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ABSTRACT 

Based on a suite of state-of-the-art high-resolution A-body simulations, we revisit the so-called 
halofit model (Smith et al. 2003) as an accurate fitting formula for the nonlinear matter power 
spectrum. While the halofit model has been frequently used as a standard cosmological tool to predict 
the nonlinear matter power spectrum in a universe dominated by cold dark matter, its precision has 
been limited by the low-resolution of A-body simulations used to determine the fitting parameters, 
suggesting the necessity of improved fitting formula at small scales for future cosmological studies. We 
run high-resolution A-body simulations for 16 cosmological models around the Wilkinson Microwave 
Anisotropy Probe (WMAP) best-fit cosmological parameters (1, 3, 5, and 7 year results), including 
dark energy models with a constant equation of state. The simulation results arc used to rc-calibratc 
the fitting parameters of the halofit model so as to reproduce small-scale power spectra of the N- 
body simulations, while keeping the precision at large scales. The revised fitting formula provides 
an accurate prediction of the nonlinear matter power spectrum in a wide range of wavenumber (fc < 
30/iMpc" 1 ) at redshifts < z < 10, with 5% precision for k < lhMpc^ 1 at < z < 10 and 10% 
for 1 < k < 10/iMpc _1 at < z < 3. We discuss the impact of the improved halofit model on weak 
lensing power spectra and correlation functions, and show that the improved model better reproduces 
ray-tracing simulation results. 

Subject headings: cosmology: theory - large-scale structure of universe - methods: N-body simulations 
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1. INTRODUCTION 

The large-scale structure of the Universe has evolved 
under the influence of cosmic expansion and gravity, and 
its statistical nature contains valuable cosmological in- 
formation. Among others, the power spectrum P(k) is 
one of the most fundamental statistical quantities char- 
acterizing the large-scale structure. It has widely been 
used for cosmological studies, both in predicting vari- 
ous observable quantities and in extractin g cosmological 
information fro m the observations (e.g.. [Peebles! [l993t 
lDodelsodl2003l K Given growing interests in high preci- 
sion cosmological observations, of particular importance 
is an accurate theoretical template of the power spec- 
trum, taking account of the nonlinear gravitational evo- 
lution. 

Weak lensing induced by the large-scale structure 
between observed galaxies and the observer provides 
a unique opportunity to directly probe matter inho 
mogeneities in the Universe. This cosmic shear sig- 
nal has been measured with a high signal-to-noise ra- 
tio by current large surveys including Canad a- France 
Hawa ii Telescope Legacy Survey fCFHTLS: iFu et al 



200 3 ), and Large Synoptic Surve ; 



LSST Science Collaborations et al 



e y lef t 



2008), 



2011 



Sloan D i gital Sky Survey (SDSS; iLin et al 
iHuff et all |20rl 



and Cosmic Evolution Sur- 



vey (COSMOS; iMassev et all 120071 : iSchrabback et~aT1 
l2Tjl0f) . These surveys provided useful constraints on 
the cosmological parameters such as the matter den- 
sity parameter Q m and the amplitude of density fluc- 
tuation er 8 . Futur e surveys such as Subaru Hyper 
Suprime-Cam C HSGlMivazaki et al.ll2006D. Dark Energy 
burvcy (DES; iThe Dark Energy Survey Collaboration! 



Telescope (LSST; 
aim at mea- 
suring the cosmic shear signal with unprecedented pre- 
cisions. While weak lensing probes matter fluctua- 
tions projected along the line-of-sight, one can ex- 
tract the redshift evolution of the fluctuations, and 

hence accurate information on dark energy, us ing a 

technique called len sing tomography (e.g., iHul Il999t 
iTakada fc J aid 120041 ) or a cross-correlation with inter- 
vening objects (e.g.. lOguri fc Takadal l20l"DY However, 
accurate and unbiased cosmological constraints from 
these lensing measurements can be obtained only if we 
have appropri ate likelihood function w ith given marginal 
distributions (|Sato et all l2010t 1201 ll) and an accurate 
model of the power sp ectrum P{k). For instance, 
iHuterer fc Takadal (|2005f ) argued that we typically need 
a few percent accuracy of P{k) at the wavenumber 
k < 10/iMpc -1 in order for the uncertainty of P(k) not 
to degra de cosmolog ical constraints in DES and LSST 
(see also I Eiflerl 1 2 1 if in which a similar conclusion is ob- 
tained). 

In the linear and quasi-linear regime of density fluc- 
tuations, the power spectrum can be computed for any 
given initial conditions and c osmological parameters us- 
ing perturbation theory fe.g-. lBernardeau et al.|[2002l for 
a review). In the nonlinear regime, however, one has 
to resort to cosmological A-body simulations to study 
the nonlinear gravitational evolution. A-body simula- 
tion results are then used to develop phcnomcnological 
halo models or fitting for mulae of nonlinear gravita tional 
clustering. For instance, iPeacock fc DoddsT (|1996l ) pro- 
vided a fitting formula of P(k) based on a scaling ansatz 
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presented in lHamilton et al.l (|1991|) . ISmith et all (|2003l 

hereafter S03) proposed a new model of P(k), the so- 
called halofit model, whic h is based on a halo model of 
struc ture formation (e.g., ISeliakl [20001 : ICoorav fc Shethl 
2002). In this halo model, all the matter content in 
the Universe is assumed to be bound in dark matter ha- 
los. Then the power spectrum is decomposed into two 
terms, the so-called one- and two-halo terms. The one- 
halo term describes matter correlations within the same 
dark matter halo, and is determined by the density pro- 
file of each halo. On the other hand, the two-halo term 
arises from the correlation between two distinct halos. 
The one-halo term dominates at small scales, whereas 
the two-halo term dominate at large scales. The halofit 
model chose the functional form of P(k) based on the 
halo model, but the model parameters were calibrated 
from TV-body simulation results. 

The halofit model by S03 is widely used to calculate 
the nonlinear matter power spectrum, yet it has been 
reported that the model fails to reproduce recent high- 
resolution TV-b o dy simulation results a t sma ll scales (e.g., 
Springel et all l2005t [H ilbert c t~aT1 120091: ISato et al] 
2009tlBovlan-Kolchin et alf2 009: Takaha shi et al. 112011b 
Kiessling et al.l 12011b IValageas fc Nishimichil l2011allbl : 
Harnois-Deraps et al.l 12012) : llnoue fc Takahashil |2012| ). 
For instance. IWhite fc Valef (|2004f ) first pointed out that 
the halofit predicts a s maller power than thei r numerical 
results at small scales. iHeitmann et al.l (|2010t) ran a suite 
of high-resolution simulations, called "Coyote Universe" , 
and showed that P(k) predicted by the halofit is ~ 5% 
smaller than their numerical results at k ~ 1/iMpc -1 . 
The reason of the difference comes from the fact that 
the TV-body simulations used in S03 have lower spatial 
resolution than latest ones. An outcome of the Coyote 
Universe simulations is a publicly available code "cosmic 
emulator" to calculate the nonlinear matter power spec- 
trum by interpolating the s imulations results for 38 dif- 
ferent cosmological models (jLawrence eEaDSOl^). How- 
ever, their emulator is restricted to a narrow range in 
k < 3/iMpc -1 and at low redshift < z < 1. Also, 
the Hubble parameter is automatically specified in the 
code using the cosmic microwave background (CMB) 
anisotropy constraint on the distance to the last scat- 
tering surface. 

In this paper, we revisit the halofit model based on 
state-of-the-art high-resolution TV-body simulations in 16 
cosmological models around the Wilkinson Microwave 
Anisotropy Probe (WMAP) best-fit cosmological param- 
eters. We allow the dark energy equation of state w to 
deviate from — 1, assuming that w does not evolve with 
redshift. While the original halofit model in S03 contains 
30 parameters, we increase the number of parameters to 
35 in order to achieve a better fit to the simulations. 
The new formula we present, which is summarized in 
Appendix, is widely applicable in the wavenumber range 
of k < 30/iMpc -1 and the redshift range of < z < 10. 
Simply replacing the parameters in the original halofit 
model with new ones in the Appendix, an accuracy of 
fitting function is improved especially at small scales. 

The present paper is organized as follows. In Section^ 
we begin by describing the TV-body simulations and cos- 
mological models used for the power spectrum analysis. 
Combining the TV-body results with different box sizes, 



TABLE 1 





n h 


S 2 m 


h 


o"8 


n s 


—10 


WMAP1 


0.044 


0.29 


0.72 


0.9 


0.99 


1 


WMAP3 


0.041 


0.238 


0.732 


0.76 


0.958 


1 


WMAP5 


0.046 


0.279 


0.701 


0.817 


0.96 


1 


WMAP7 


0.046 


0.272 


0.7 


0.81 


0.97 


1 


WMAP 7a 


0.046 


0.272 


0.7 


0.81 


0.97 


0.8 


WMAP7b 


0.046 


0.272 


0.7 


0.81 


0.97 


1.2 



Note. — Best-fit cosmological parameters in a series of 
WMAP papers. Here we show the baryon density Q^, the 
matter density f2 m , the Hubble constant h, the amplitude 
of power spectrum at Sh~ 1 Mpc erg, the spectral index n s , 
and the equation of state of dark energy w. We assume a 
flat curvature (f2 w = 1 — C m ). 



TABLE 2 





n h 


S2m 


h 


C8 


n B 


— w 


mOO 


0.0432 


0.25 


0.72 


0.8 


0.97 


1 


mOl 


0.0647 


0.4307 


0.5977 


0.8161 


0.9468 


0.816 


m02 


0.0637 


0.4095 


0.5907 


0.8548 


0.8952 


0.758 


m03 


0.0514 


0.2895 


0.6763 


0.8484 


0.9984 


0.874 


m04 


0.0437 


0.2660 


0.7204 


0.7 


0.9339 


1.087 


m05 


0.0367 


0.2309 


0.7669 


0.8226 


0.9726 


1.242 


m06 


0.0462 


0.3059 


0.7040 


0.6705 


0.9145 


1.223 


m07 


0.0582 


0.3310 


0.6189 


0.7474 


0.921 


0.7 


m08 


0.0428 


0.2780 


0.7218 


0.8090 


0.9855 


1.203 


m09 


0.0623 


0.3707 


0.6127 


0.6692 


0.979 


0.739 



Note. — Cosmological parameters of Coyote models. 



we discuss in detail the convergence of power spectrum 
measurement over the wide range of wave number. In 
Section [3J we re-calibrate the halofit model, and the re- 
vised version of the halofit model, whose explicit formula 
is given in Appendix, is compared with our TV-body sim- 
ulations. As an important implication of revised halofit 
model, in Section [4j we compute weak lensing power 
spectra, and compare them with direct ray-tracing sim- 
ulation results, particularly focusing on the small-scale 
behavior. Finally, Section [S] is devoted to conclusion and 
discussion. 

2. N-BODY SIMULATIONS 

2.1. Power Spectrum 

In this section, we describe our cosmological TV-body 
simulations used in this paper. We follow the nonlin- 
ear gravitational evolution of 1024 3 collisionless particles 
in a cubic box of side L. We use the public cosmolog- 
ical TV-bo dy simulation code Gadget2 which is a tree- 
PM code ([Springe! et al.l 120011: ISpringell l2005ft . We use 
2048 3 PM grid to follow the gravitational evolution at 
small scales accurately. We generate the initial condi- 
tions based on the second-order Lagrangian perturba- 
tion theory f2LPT: IcTocce et al.l [20061: iNishimichi et al.l 
l2009h with the initial linear power spectrum calculated 
by the Code for Anisotropics in the Microwave Back- 
ground (CAMB: iLewis et al.ll2000h . The initial redshift 
is set to z- m = 99. We store simulation results (particle 
positions) at various redshifts from z = to 10. The 
softening length is set to 5% of the mean particle sepa- 
ration. To calculate the power spectrum, we assign the 
particles on TV 3 = 1280 3 grid points using the cloud- 
in-cells (CIC) method (jHocknev fc Eastwoodl 119811 ) to 
obtain the density field. After performing the Fourier 
transform, we correct the window function of CIC by di- 
viding each mode by the Fourier transform of the window 
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TABLE 3 

L(fe~ 1 Mpc) JVg Nr fcNyq^Mpc' 1 ) r^fe-^pc) z init z ou t 
WMAP models 2000 1024 3 3 O 97S 99 0, 0.35, 0.7, 1, 1.5, 2.2, 3, 5, 7, 10 

800 1024 3 3 4.0 39.0 99 0, 0.35, 0.7, 1, 1.5, 2.2, 3, 5, 7, 10 

320 1024 3 6 10. 15J5 99 0, 0.35, 0.7, 1, 1.5, 2.2, 3, 5, 7, 10 

Coyote models 1000 1024 3 1 3l 50 99 0, 0.35, 0.7, 1, 1.5, 2.2, 3, 5, 7, 10 

320 1024 3 1 10. 15S 99 0, 0.35, 0.7, 1, 1.5, 2.2, 3, 5, 7, 10 

Note. — Model parameters of our numerical simulations for the WMAP models (upper rows) and the Coyote models 
(lower rows): the box size L, the number of particles JV 3 , the number of realizations JV r , the Nyquist frequency &Nyq = 
(2-7T / L)(N p /2), the softening length r B , the initial redshift Zi n i t and the redshifts of the simulation outputs z ou t. 



kernel as S k -> Y\i=x,y,z [sinc(I/fci/2iVg)p 2 x <S k , where (5k 
is the density fluctuation in Fourier space and sinc(x) = 
sin(x )/x (e.g.. iTakahashi et al.l l2009: Sat o fc Matsubaral 
1201 ll ). In addition, to evaluate the power spectrum at 
small scales accurately, we fold the particle positions 
into a smaller box by replacing x — > x%(L/2 n ) where 
the operation a%b stands for the remainder of th e divi- 
sion of a by & (e.g., IValageas fc Nishimicm]|2011al ). This 
procedure leads to effectively 2™ times higher resolution. 
Here we adopt n = 0, 2, and 4. We use the density 
fluctuation 5k up to half the Nyquist frequency deter- 
mined by the box size L/2 n with the grid number N g , 
i.e., 1/2 x fc Nyq = (ir/(L/2 n ))(N g /2), with n = 0, 2, 
and 4. This condition corresponds to fc < 6.3, 25, and 
100/iMpc -1 with n = 0, 2, and 4, respectively, for the 
box size of L = 320/j -1 Mpc with N g = 1280. Finally, we 
compute the power spectrum 

iv fc k 

where the summation over Fourier modes is done for the 
modes falling into the bin [fc — Afc/2, fc + Afc/2], and 
Nk denotes the number of available Fourier modes in the 
bin. We do not subtract the shot noise in the measured 
power spectrum. Instead, we do not use P(fc) at small 
scales where the shot noise dominates (see Section [3]) . 

2.2. Cosmological Models 

In this paper, we use simulation results for 16 cosmo- 
logical models. Six are taken from the results of WMAP 
papers and 10 are from the cosmological models adopted 
by the Coyote Universe. For all of the models, we assume 
a flat curvature (fi w = 1 — f2 m , where fl w is the dark 
energy density). The first four WMAP models, which 
are shown in Table [TJ are the bes t-fit ACDM models 
of WMAP 1, 3, 5, 7 year resu lts (|Spergel et al.l I2003L 
l2007tlKomatsu et a l. 2009, 20lTl). The other two models, 
WMAP 7a and WMAP 7b, are the same as the WMAP 7 
model except that we slightly change the equation of 
state parameter of dark energy (w = —0.8 and —1.2). 

In addition, we also examine 10 models among 38 cos- 
mological models presented in the Coyote Universe, as 
shown in Table [5J These models, tagged as mOO to m09, 
were use d in a series of papers of the Coyote Universe 
proje ct (jHeitmann et al.1 I2009L 120101 : lLawrence et al.1 
|2010j), in which 38 cosmological models in a parame- 
ter range of 0.120 < Q^h 2 < 0.155, 0.0215 < Q h h 2 < 
0.0235, 0.85 < n s < 1.05, -1.30 < w < -0.70, and 
0.61 < erg < 0.9 are used to make a fitting function 
of the power spectrum. The models used in our paper 
correspond to their first 10 cosmological models. Since 



the cosmological parameters in Coyote models arc differ- 
ent from the WMAP models typically by 20 - 30%, we 
use our simulation results for these Coyote cosmological 
models to check the dependence of our fitting function 
on cosmological parameters. 

Table summarizes our simulation setting, including 
the box size, the number of particles, the number of 
realizations, and the softening length, for the WMAP 
and Coyote models. In the WMAP models, we adopt 
the simulation boxes of L = 2000, 800, and 320/i _1 Mpc 
on a side. We prepare 3 different random realizations 
for L = 2000 and 800/i~ 1 Mpc, and 6 realizations for 
320ft~ 1 Mpc to reduce the sample variance. We combine 
the power spectrum P(k) obtained from the different 
simulation boxes to cover a wide wavenumber range. The 
specific procedure for combining the -P(fc) is discussed in 
the next subsection. We use the mean power spectrum of 
these realizations. The Nyquist wavenumber of the mean 
particle separation in the smallest box (L = 320/i _1 Mpc) 
is fc = 10/iMpc -1 . The bin width is set linearly, Afc = 
0.01/iMpc _1 , in the linear regime fc < 0.3/iMpc -1 , and 
logarithmically, Alog 10 (fc//iMpc _1 ) = 0.02, in the non- 
linear regime fc > 0.3ft. Mpc" 1 . We analyze ten outputs 
at redshifts < z < 10, z = 0, 0.35, 0.7, 1, 1.5, 2.2, 3, 5, 
7, and 10. We checked that the power spectra of our 
simulations agree with the results of higher resolution 
simulations, in which we set the finer simulation param- 
eters for the time step, the force calculation, etc., within 
2(6)% for fc < 10(30)ftMpc" 1 . 

In the Coyote models, we use the simulation boxes 
of lOOO/i^Mpc and 320/i _1 Mpc on a side. We pre- 
pare a single realization for each of the 10 models. 
The Nyquist wavenumber is fc = 3.2(10)/iMpc _1 for 
L = 1000(320)/i _1 Mpc and we use a logarithmic bin of 
Alog 10 (fc/ftMpc _1 ) = 0.1. We use ten outputs at the 
redshifts < z < 10, z = 0, 0.35, 0.7, 1, 1.5, 2.2, 3, 
5, 7, and 10. We confirmed that our simulation results 
agree with the cosmic emulator results within 3% for 
0.1 < fc < 3/iMpc -1 at < z < 1. 

Finally, as a further cross check, we compare our simu- 
la tion results with two high-resolu tion simulation results 
inlValageas fc Nishirmchir()2011a| ) and ITakahashi et al. I 
(|2011f ). In both the simulations, the same codes, 
Gadget-2 with 2LPT initia l condition, were used. 
IValageas fc Nishimichil (|2011af ) employed 2048 3 particles 
in different box sizes of L = 4096, 2048, 1024, and 
512/i _1 Mpc. They calculated the power spectra at red- 
shifts z = 0.35, 1, and 3 in the WMAP5 model. The 
Nyquist wave numbe r in the smallest box is l S/tMpc -1 . 
On the other hand, ITakahashi et al. 1 (|2011h employed 
1024 3 particles on the box size of 50/i _1 Mpc at redshift 
z = — 20. They prepared independent four realizations. 
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Fig. 1. — Power spectra P(k) from the different simulation 
box sizes for the WMAP5 model at redshifts z = 0.35, 1, 3, 5, 
and 10. The vertical axis shows power spectra in our simula- 
tions normalized by the theoretical nonlinear matter power spec- 
trum model from the original halofit model in S03. Green, blue, 
and orange symbols are simulation results from the box sizes of 
L = 2000, 800, and 320/i _1 Mpc, respectively. Red symbols arc 
the same as the oranges, but including the folding method with 
n = 2 (see S ection I2.lt. Gray sy mbols show simulation results 
from lValageas fe~ Nishiinichi (2011a). Arrows denote wavenumbers 
where P(k) with the different box sizes are connected. Dashed 
curves denote the shot noise. Vertical solid lines indicate the max- 
imum wavenumber fc max for deriving our fitting formula. 



The cosmological parameters are based on the WMAP 
5 year result, although the values were slightly different 
from those of the WMAP5 model listed in Table 1. We 
use 14 outputs at redshifts z = — 10, z = 0, 0.35, 
0.7, 1, 1.5, 2.2, 3, 3.4, 4.1, 4.7, 5.5, 6.4, 7.6, and 9.1. 
The Nyquist wave number is 64ft Mpc" 1 . We use these 
simulation results in our analysis for only small scales, 
fc > lO^ftMpc" 1 at z < 3(> 3), in order to check the 
asymptotic behavior of our fitting formula at high fc limit. 

2.3. Accuracy of Our N-body Simulations 

In this subsection, we compare our simulation results 
with previous works to check the accuracy and conver- 
gence of our A^-body simulations. We also describe the 
procedure for combining the P(k) from the different sim- 
ulation box sizes. Figure [T] shows the power spectrum 
P(k) of our simulations for the different box sizes for 
the WMAP5 model at z = 0.35, 1, 3, 5, and 10. The 
vertical axis shows the measured power spectra normal- 
ized by the theoretical model of the nonlinear power 
spectra from the original halofit model in S03. Green, 
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Fig. 2. — Power spectra P(k) from our simulations of L = 
1000/i~ 1 Mpc divided by those of the Cosmic emulator for the ten 
Coyote models at z = (top) and 1 (bottom). Red symbols with 
error bars are for the fiducial model mOO. Horizontal dotted lines 
indicate the fractional error of 3%. 



TABLE 4 



z 


fc 20(>(>/8()() 


^800/320 


^min ^max 


0,0.35 


0.1 


1. 


0.02 30 


0.7,1 


0.1 


0.6 


0.02 20 


1.5,2.2 


0.1 


0.1 


0.02 10 


3 


0.1 


0.3 


0.02 8 


5 


0.1 


0.3 


0.02 6 


7 


0.1 


0.3 


0.02 3 


10 


0.1 


0.3 


0.02 2 


Note. — The summary of wavenumbers (in 
units of h Mpc ) where the simulation results 



P(k) with the different box sizes are connected 
at each redshift, ^2000/800 * s f° r L = 2000 and 
800/i Mpc, and fcgnn/320 ' s f° r L = 800 and 
320/i — 1 Mpc. The minimum and maximum 
wavenumbers used in our analysis are shown 



by fcmin and 



respectively. 



blue, and orange symbols are the results from the dif- 
ferent box sizes of L = 2000, 800, and 320/i" 1 Mpc, re- 
spectively. Red symbols are the same simulation results 
as plotted in the orange symbols (L = 320ft~ 1 Mpc), 
but using the folding method with n = 2 in comput- 
ing P{k) (see Section 12. ip . We plot the mean P(k) 
with the error bars among the realizations. Gray sym- 
bols at z = 0, 0.35, and 1 are the simulation results of 
IValageas fc Nishimichil (|2011al) . Vertical arrows indicate 
the wavenumbers at which we connected the simulation 
results from the different box sizes. For example, we 
used the results from L = 800ft _1 Mpc between blue and 
the orange arrows. We connect the results of L = 800 
and 2000ft" 1 Mpc at fc = 0.1ft Mpc" 1 for all the red- 
shifts. This is because at the baryon acoustic oscillation 
(BAO) scale (0.1ft, Mpc -1 < k < 0.3ft, Mpc -1 ), the simu- 
lation results with L = 800/i _1 Mpc show a better agree- 
ment with tho se obtained from the i mprov ed perturba- 
tion theory by iTaruva fc Hiramatsul (|2008f ). The initial 
redshift of z init = 99 is high for L = 2000ft" ^pc with 
1024 3 particles, and hence P{k) from L = 2000ft~ 1 Mpc 
at low redshifts is slightly smaller than the -P(fc) from 
L = 800ft" 1 Mpc. Next, in connecting the P{k) from 
L = 800ft" 1 to 320ft," 1 Mpc, we use P{k) from L = 
800ft" : Mpc up to the wavenumber where P(k) agrees 
with the previous high-resolution simulation results (gray 
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symbols) at z = 0.35, 1, 3. In this way, the connecting 
scales arc determined to k — 1, 0.6, and 0.3ft Mpc" 1 for 
z = 0.35, 1, and 3, respectively. For the redshifts z = 0, 
0.7, 1.5, and 2.2, we interpolate the scales derived above, 
and for the redshifts z = 5, 7, and 10 we simply adopt 
the same result as that at z — 3. The connecting scales 
at each redshift are summarized in Table |U While these 
connecting scales are derived from the WMAP5 model, 
we also use the same connecting scales in Table 0] for the 
other cosmological models. We confirmed that the power 
spectra from different simulations connect smoothly at 
these scales in all the cosmological models studied in this 
paper. 

Figure [2] shows the power spectrum P(k) in our simula- 
tions of L = lOOOftMpc -1 divided by that of the Cosmic 
emulator for the ten Coyote models at z = 0, 1. The col- 
ored symbols correspond to the ten Coyote models. Red 
symbols with error bars are for the fiducial model mOO. 
Here the error bars show the Gaussian errors because 
we have only one realization for each Coyote cosmologi- 
cal model. As clearly seen in the Figure, our simulation 
results agree with the Cosmic emulator within 3% for 
fc = 0.1- 2ft Mpc" 1 . 

3. HALOFIT MODEL 

In the halofit model, the power spectrum consists of 
two terms (S03): 



A 2 (fc) 



A H (fc), 



(2) 



where A 2 (fc) = fc 3 P(fc)/(27r 2 ) is the dimcnsionless power 
spectrum. The first term is called the two-halo term 
that dominates at large scales, whereas the second term 
is referred to as the one-halo term that is important at 
small scales. We adopt almost the same functional form 
as in S03 for both the two terms in Equation ([2]) . We use 
our high-resolution simulation results to re-calibrate the 
model parameters of the halofit formula so as to minimize 
the discrepancies. To do so, we employ the standard chi- 
squared method to find the best-fit solution: 



X 



E E E^ 



[-Pi, model z ) Pj,sim {k, z)J 



i k=k m in z — 

(3) 

where P^modci is the model prediction, Pi, s im is the sim- 
ulation results, and i runs over the WMAP cosmological 
models shown in Table [T] We note that we use only 
the six WMAP models in this chi-squared analysis. The 
remaining ten Coyote models are used to check the accu- 
racy of our fitting formula. We simply set the variance 
of = -Pj 2 s j m , and the weight function is set as follows: 

W(k, z) = 10, k < 0.3/i Mpc" 1 & < z < 3, 

= 1, 0.3/iMpc" 1 < k < 10ft Mpc" 1 & < z < 3, 
= 0.2, k < lOftMpc" 1 & 3 < z < 10, 
= 0.1, k > 10ft Mpc -1 . 

The weight factor is chosen so that the final fitting for- 
mula gives a better accuracy at the BAO scales at low 
redshifts. 

In (quasi-)linear regime, the error of P(k) is given by 
the Gaussian error which is th e inverse of the square root 
of the number of modes (e.g.. iFeldman et al.lll994| ). We 



consider only the wavenumber bins where the Gaussian 
error of P(k) is less than 3% for the fitting, which corre- 
spond to fc m in = 0.02ft Mpc" 1 . In nonlinear regime, the 
non-Gaussian error arises due to the mode coupling, bu t 
it is smaller than 5% (see e.g., Takaha shi et al.l I2009D . 
While the Nyquist wavenumber is k = 10ft Mpc" 1 , we 
sum up the wavenumber up to k = 30ft Mpc -1 , be- 
cause the nonlinear power spectrum is reliable down 
to scales correspond ing to the softening length (e.g., 
lHamana et al.l [2002t) . On the other hand, we do not 
use the wavenumber where the shot noise dominates 
the power spectrum. The maximum wave number is 



fcmax = 30ft Mpc" 
at z = 0.7 and 1, k„ 



and kn 



at z — 0, 0.35, fc max 
= 10ft Mpc" 1 at z 



8, 6, 3, and 2ft Mpc" 1 at 



= 20ftMpc -1 
1.5 and 2.2, 
3, 5, 7, and 



10, respectively. The minimum and maximum wavenum- 
bers (fc m in and fc max ) are listed in Table 2) For all the 
WMAP models, the power spectrum at k = k max is 10 
(3) times larger than the shot noise at z < 3 (5 < z < 10). 
There are 35 free parameters in our revised halo-model 
(30 parameters in the original model). We summarize 
the best-fit parameters in Appendix. 

Figure [3] shows the power spectra P(k) as a function 
of the wavenumber k for the WMAP1, 3, 5, and 7 mod- 
els at redshifts z = 0, 0.35, 1, and 3. Here, to empha- 
size the difference between the simulation results and the 
theoretical models, the power spectrum P(k) is multi- 
plied by the factor k 15 . Black circles with error bars are 
our simulation results, whereas gr ay symbols in WMAP5 
are the simulation results from IValageas fc Nishimichil 
(I2011al) . As seen in the Figure, our simulation results 
agree with the results in IValageas fc Nishimichil (|2011al) 
very well. Red curves show our fitting function, which are 
significantly better than the original halofit model in S03 
shown by black curves. As clearly seen in the Figure, the 
original halofit model grossly underestimates the power 
spectra at smaller scales (k > 0.1ft Mpc -1 ). Our model 
agrees with the simulation results very well down to small 
scales for all the cosmological models. The agreement of 
our fitting formula with simulations is better than 8% at 
k < 10ft Mpc" 1 for the WMAP cosmological models in 
the redshift range of < z < 3. For all the WMAP cos- 
mological models at k m - m < k < /c max , the rms deviation 
of our best-fit model from the simulation results is 1.7% 
at < z < 10. 

Figure [4] shows the same results as in Figure |3l but we 
focus on the results at the BAO scale of 0.05ft Mpc" 1 < 
k < 0.3ft Mpc" 1 . In the vertical axis, the power 
spectrum is normalized by the smooth nonlinear power 
spectrum P^ 3 (k), which is calculated by using a no- 
wiggle fitting formula of lEisenstein fc Hul (|1998f ) with 
nonlinear corrections computed by the original halofit 
model in S03. Green curves show theoretical predic- 
tions obtained from the improved perturbation theory 
called closure theory, which efficiently resumms a class 
of infi nite series of higher-order perturbative correc- 
tions dTaruva fc Hiramatsu 120081 : iHiramatsu fc Taruval 
120091 : iTaruva et al.ll2009[ ). These predictions include the 
corrections at the 2-loop order based on the Born approx- 
imation. The Figure indicates that our model agrees with 
the simulation results better than S03 especially at low 
redshifts. However, the closure theory shows even bet- 
ter agreements in the quasi-linear regime. In the BAO 
scales, our fitting formula reproduces simulation results 
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10" 1 10° 10 1 10" 1 10° 10 1 10 1 10° 10 1 10 1 10° 10 1 



k(h/Mpc) 

Fig. 3. — Power spectra P(k) for the WMAP cosmological models at z = 0, 0.35, 1, and 3. In the vertical axis, the power spectrum 
P(k) is multiplied by the factor k in order to show the differences between the simulation results an d the theoretical model s clearly . 
Black filled circles with the error bars plot our simulation results, and gray symbols are the results from Valageas & Nishimichi (2011a). 
Red solid curves show our revised halofit model (see Appendix), whereas black solid curves show the original halofit model in S03. Black 
dashed curves plot the linear power spectra. 



within 4.6% for the WMAP models at < z < 3. 

Figure [5] shows ratios of the measured power spec- 
tra in our simulations to the revised fitting formula, 
-P s im(fc)/-Pmodci(fc), for all the cosmological models at 
z = 0, 1, 3, 5, and 10. Filled circles shows results for 
the six WMAP models, whereas plus and cross sym- 
bols are for the ten Coyote models. The plus and cross 
symbols are the results from larger (smaller) box sizes 
L = 1000(320)/T. _1 Mpc and are shown only for larger 
(smaller) scales of k < 1(> 0.3)ft-Mpc _1 . The horizontal 
dotted lines show the errors of 5%. Vertical solid lines at 
z = 1, 3, 5, and 10 indicate the maximum wavenumbcr 
fc max in the chi-squared calculation in Equation ([3]). Note 
that fcmax = 30ft. Mpc" 1 at z = 0. As seen in the Figure, 



relative errors are typically less than 5(10)% at z < 3(10) 
for all the cosmological models. Although we did not in- 
clude the Coyote models in our fitting, our fitting formula 
reproduces simulation results for the Coyote models very 
well, with the errors less than 8% for k < lh Mpc -1 at 
z = - 10, and 13% for k < lO/iMpc" 1 at z = - 3. 
Our fiducial models (WMAP5 and 7) show better agree- 
ment, while the other models involves slightly larger er- 
rors. At z = 0, the cosmological models with dark en- 
ergy (w ^ —1) and with high erg (WMAP1) show larger 
errors. The cosmological models with large (small) equa- 
tion of state w show the larger (smaller) simulation re- 
sults Psim(k) than our fitting model. For example, for 
the Coyote mOl (w = -0.816) and m02 (w = -0.758) 
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WMAP1 WMAP3 WMAP5 WMAP7 




j i i i i i i i ' + i i i i i i i i i + i i i i i i i i i + i i i i i i i i i_j 

0.1 0.2 0.3 0.1 0.2 0.3 0.1 0.2 0.3 0.1 0.2 0.3 



k(h/Mpc) 

Fig. 4. — Similar to Figure^ but focusing on the BAO scale of 0.05ft Mpc -1 < k < 0.3ft Mpc^ 1 . The vertical axis show P(k) normalized 
by the smooth nonlinear power spectrum P 1 ^ 3 (fc) which is calculated by using a no- wiggle fitting formula of Eiscnstcin & Hu (f998) and 
the original halofit model in S03 for the nonlinear correction. Again, red solid curves show our revised halofit model, whereas black solid 
curves sh ow the original halofit model in S03. Green curves show the prediction by the closure theory which is a higher-order perturbation 
theory by Taruva & Hiramatsu (2008). Dashed curves show the linear power spectra. 



models, shown as the brown and orange crosses, the sim- 
ulation results are over 10% larger than our best fitting 
model for k ~ lOft-Mpc -1 at z = 0. At higher redshifts 
z > 3, the increase of the ratio at high k are due to the 
shot noise. At high z, the errors depend mainly on the 
spectral index n s since cosmological models converge to 
the Einstein de-Sitter model. The models with the steep 
(shallow) spectral index shows the large (small) ratio at 
small scales k > 1/iMpc^ 1 . 

Peacock also provided an improved halofit model which 
gives simply factor two times larger power than the 
original model for small scales k > 10/iMpc -1 (see his 



homepage 1 ). But, his model predicts a smaler power 
than our simulation results at quasi-linear scale (k = 
0.1 — 1/;, Mpc -1 ). As clearly seen in the figure [TJ the ra- 
tio of the simulation results to the halofit in S03 is not 
two for k > 10/iMpc -1 . Rather, the ratio is functions of 
the redshift, the wavenumber and the cosmological pa- 
rameters. 

4. IMPLICATIONS FOR WEAK LENSING PREDICTIONS 

In this section, we study how the revised model of the 
matter nonlinear power spectrum affects weak lensing 
observables. Specifically, we calculate the convergence 



1 http:/ /www. roe. ac.uk/~jap/haloes/ 
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Fig. 5. — Our simulation results of the power spectrum P s i m (fc) 
divided by our improved fitting formula P mo d e i(fc) (see Appendix) 
for all the cosmological models at z = 0, 1, 3, 5, and 10 (from top 
to bottom). Filled circles show results for the six WMAP models, 
whereas plus and cross symbols are for the ten Coyote models. 
Blue and red big circles are WMAP5 and WMAP7. Black, green, 
gray and yellow circles arc results for the WMAP1, 3, 7a and 7b 
models, respectively. Horizontal dotted lines indicates the relative 
error of 5%. 

power spectra, correlation functions, and the CMB tens- 
ing using the new fitting formula that is summarized in 
Appendix. We compare results based on our fitting for- 
mula with those from the original halofit model in S03 
as well as the direct ray-tracing simulation results. 

Images of distant galaxies are distorted by grav- 
itational l ensing due to intervening matter fluctua- 
tions (e.g.. iBartelmann fc Schneider! 1200 It iMunshi et al.l 
120081 for a review). The image deformation is char- 
acterized by the lensing convergence K and shear 7. 
The convergence field is expressed as the integration of 
a weighted three-dimensional density fluctuations along 
the line-of-sight. Hence, the convergence power spectrum 
Ce can be expressed as a projection of the matter power 
spectrum weighted with the radial lensing kernel along 
the line-of-sight: 

C l = J & dzW ( z ' z ^ P { k= ^j' z )> (4) 

where W(z, z s ) is a weight function. 

Figure |5] shows the convergence power spectra and 
correlation functions in the WMAP3 model at source 
rcdshifts z s = 0.8 and 1.5. In both panels, red solid 
curves show the prediction using our revised model 
of the power spectrum, and black solid curves show 



that from the original halofit model in S03. Filled 
circles with error bars plo t direct ra y -traci n g sim - 
ulation results obtained in ISato et al.l (|2009l . 1201 lh . 
They used a standard ray-tracing me thod using code 
"RAYTRIX" (jHamana fc Mellieill200l . With the 256 3 
particles in the rectangular box of 240 and 480/i _1 Mpc 
on each side, they prepared 1000 convergence and shear 
maps in the field of view 5x5 deg 2 . The mean and 
error shown in Figure [6] are estimated from the 1000 re- 
alizations. Note that the size of error bars is inversely 
proportional to the square root of the sky coverage. 
For example, the Subaru HSC wide survey will observe 
~ 1500 deg 2 , which suggests that expected error bars 
are ~ 8 times smaller than those plotted here. Vertical 
arrows indicate the multipole below which their simu- 
lations are consistent with higher resolution simulations 
(512 3 particles) within 5%. The Figure indicates that our 
model predictions agree with the simulation results much 
better than those of the original halofit model. This 
suggests that the use of the improved fitting formula as 
presented in this paper is essential to extract cosmologi- 
cal information from future high-resolution weak lensing 
measurements . 

We also consider the weak lensing effect on the CMB 
temperature anisotropy. Gravitational lensing by fore- 
ground matter distributions is known to affect the light 
path of CM B photons coming from the last scattering 
surface (e.g.. iLewis fc Challinorl [20061 ) . As a result, the 
spatial pattern of CMB temperature fluctuations is dis- 
torted, which leads to the modification of the tempera- 
ture power spectrum. Since the lensing deflection angle 
is typically a few arcmin, the temperature power spec- 
trum shape is significantly modified at smaller scales of 
I > 1000. Figure [7] shows the predicted lensed CMB 
temperature power spectrum. We use the CAMB code 
to compute the unlensed (primordial) CMB power spec- 
trum, which is plotted by the dashed curve. The red 
solid curve is the lensed power spectrum using our re- 
vised fitting formula, while the black solid curve is the 
one derived from the original halofit model in S03. The 
fitting formula of the matter power spectrum is used to 
calculate the power spectrum of the deflection angle by 
the foreground matter distribution. We find that our 
model enhances the power at small scales I > 4000 by 
- 10%. 

5. CONCLUSION AND DISCUSSION 

The halofit model presented in S03 has widely been 
used as a standard cosmological tool to predict the non- 
linear matter power spectrum. However, it has been ar- 
gued that the halofit model fails to reproduce recent high- 
resolution simulation results such that it underestimates 
the power spectrum by a few ten percent at small scales 
(fc > 0.1/iMpc _1 ). The difference is crucial for analysis 
of upcoming weak lensing surveys such as Subaru HSC 
survey, DES, and LSST. In this paper, we have revisited 
the halofit model using the high resolution simulations 
for the 16 cosmological models around the WMAP best- 
fit cosmological parameters, including the variation in 
the equation of state of dark energy. The revised fitting 
formula can reproduce the simulation results very well 
in the range of fc < 30ft. Mpc -1 and < z < 10. Our 
new fitting formula is summarized in Appendix, which 
can easily be updated from the original halofit model by 
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I 9 (arcmin) 

Fig. 6. — Convergence power spectra (left) and correlation functions (right) at the source redshifts of z s = 0.8 and 1.5. Red curves 
show theoretical predictions using our revised halofit model, and black curves s how t hose from the original halofit model in S03. Filled 
circles with error bars plot ray-tracing simulation results from [Sato ct al , (2009, 2011). Vertical arrows show the scale down to which the 
ray-tracing simulation results are accurate within 5%. Error bars are for a 5 X 5 deg 2 survey, and scale as the inverse square root of the 
ficld-of-view. 
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Fig. 7. — The lensed CMB temperature power spectrum. The 
red curve shows the prediction of the our revised halofit model, 
while the black curve is that of the original halofit model in S03. 
The dashed curve shows the unlensed temperature power spectrum. 
Note that we plot the temperature power spectrum multiplied by 
£ 6 , i.e., l 6 C t /(2ir), not the usual convention of t(t + l)C/>/(27r). 



simply replacing parameters in original model with new 
values as well as adding a few terms. 

We comment on effects of baryon cooling and massive 
neutrinos, both of which affect P(k) at small scales. The 
baryon cooling would enhance the power at small scales 
by some ten percent at k = 10/iMpc _1 and the enhance- 
ment becomes more significant fo r smaller scales (e.g. , 
Jing et al.l 120061: [Rudd et~aTl l2008l ICasarini et all 120111: 
van Daalen et al. 1 2 1 If) . However, the reliability of sim- 



ulation results strongly relies o n galaxy formation mod - 
els they adopted. For example, Ivan Daalen et all (|2011l ) 
showed that the AGN feedback can decrease the power 
spectrum by > 10% at k > 1/iMpc -1 . The mas- 
sive neutrinos also suppress the growth of density fluc- 
tuation below the so-called fr eestreaming scales (e.g., 
iBrandbvee k Hannestadl l2009t iBird et all 12012ft . The 



power spectrum is suppressed by a few ten percent at 
small scales of k ~ 0.1/iMpc -1 , depending on the total 
mass of neutrinos. Even though these effects can modify 
the small-scale nonlinear matter power spectrum, an ac- 
curate knowledge of the original (dark matter only) non- 
linear power spectra as presented in this paper is still 
important as an ingredient for building models of more 
realistic nonlinear power spectra which take these effects 
into account. 

Finally, while we have improved the fitting formula us- 
ing the simulations, there are also several attempts to im- 
prove the halo model analytically. For instance, combin- 
ing the perturbation t heory at large scales wit h the ha lo 
mo del at small sca l es, IValageas k Nishimichil (|2011al Tbl) 
and lValageas et all (|2012alp]) presented an imp roved halo 
model. On smaller scales. iGiocoli et al.l (12010D provided 
a prediction of the power spectrum using the halo model 
including the effect of substructure in the individual halo. 
These models also reproduce the simulation results well 
and are consistent with our fitting formula. 
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APPENDIX 

FUNCTIONAL FORM OF THE REVISED HALOFIT MODEL 

In this Appendix, we provide the functional form of the revised halofit model. The nonlinear power spectrum, 
A 2 (fc) = k 3 P(k)/(2ir 2 ), consist of one- and two-halo terms: 

A 2 (k) = A 2 Q (k) + A 2 i (k). (Al) 

The two-halo term Aq(/c) is given by, 



A 2 (fc) = A 2 (fc) 



|i + A 2 (fc)} A ' 

l + a„A 2 (fc) 



e - f{v \ (A2) 



where A£ = k Pi,(k)/(2ir 2 ), f(y) = y/A + y /8, and Ph(k) is the linear power spectrum. The one-halo term Afj(A:) is 
written as 

A|(fc) = 7 with A' 2 (fc) = ™ = , (A3) 

where y is the dimensionless wavenumber, y = k/k a . The nonlinear scale k~ l is defined by 

a 2 (k- 1 ) = l with a 2 (R) = J d\nk Al{k)e~ k2R2 . (A4) 
The effective spectral index n e ff and the curvature C are defined as 



"off + 3 = - 



d In er 2 (i?) 



rflni? 



C = - 



d 2 \na 2 {R) 



dlni? 2 



(A5) 



The parameters a n , b n , c n , 7 n , a n , /3 n , /i n , and v n in Equations (|A2[) and (|A3|) are given by polynomials as a functions 
of n e ff and C . We determine the coefficients in the polynomials by fitting the model to our simulation results, as 
described in Section [3] The best- fit parameters are 

log 10 a n = 1.5222 + 2.8553?i ofr + 2.37067i 2 ff + 0.9903r^ ff + 0.2250n^ ff - 0.6038C + 0.1749ft w (z) {1 + w) , (A6) 

log 10 b n = -0.5642 + 0.5864n off + 0.5716n 2 ff - 1.5474C + 0.2279ft w (z) [1+w), (A7) 

log 10 c n = 0.3698 + 2.04047i off + 0.816l7i 2 ff + 0.5869C, (A8) 

7ll = 0.1971 - 0.084377. of f + 0.8460C, (A9) 

a n = |6.0835 + 1.3373n off - 0.1959n 2 ff - 5.5274C| , (A10) 

Pn = 2.0379 - 0.7354n off + 0.3157n 2 ff + 1.2490n^ + 0.3980< ff - 0.1682C, (All) 

/in = 0, (A12) 

logio y n = 5.2105 + 3.6902n off , (A13) 



where f2 w (z) is the dark energy density parameter at redshift z. The last terms in Equations (|A6 ) and (|A7[) represent 
small correction terms for dark energy 10/-I. We use the absolute value of a n in Equation (|A10[) to avoid divergence 



in the two-halo term (Equation (j A2[) ) . Finally, /i,2.3(^) m Equation (|A3p are the same as in S03 

/l(fim) = tt ra °-° 307 , AW = ^ m °-° 585 , fs(n m ) = ^r 743 , (A14) 

where Q m is the matter density parameter at redshift z. 

In summary, one can easily revise the original halofit model in S03 by replacing the parameters a n , b n , c n , 7 n , ct n , 
P n , /i n , and v n to those listed in Equations (|A6|) - (|A13|) . 
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